package guassianElimination;

import java.text.DecimalFormat;
import java.util.Arrays;
import java.util.List;

import Jama.Matrix;

public class LuDecompostion {
	private double[][] LU;
	private int m, n, pivit;
	private int[] pivitrow;
	
	public Matrix guassianElilmination(Matrix a){
		
		LU = a.getArrayCopy();
		m =a.getRowDimension();
		n = a.getColumnDimension();
		pivitrow = new int[m];
		for(int i = 0; i< m;i++){
			pivitrow[i]=i;

		}
		pivit =1;
		double[] LUrowi;
		double[] LUcolj = new double[m];
		
		for(int j =0;j<n;j++){
		
			//	Make a copy of the j-th column to localize references.	
			for(int i=0; i<m;i++){
				LUcolj[i]=LU[i][j];
			}
			
			for(int i=0; i<m;i++){
				LUrowi=LU[i];
				
				int kmin=Math.min(i, j);
				double s=0.0;
				for(int k=0; k< kmin;k++){
					s =s + LUrowi[k]*LUcolj[k];
				}
				
				LUcolj[i] = LUcolj[i]-s;
				LUrowi[j]=LUcolj[i];
				
			}
			/**
			 * find the pivit and swap
			 */
			int p =j;
			for(int i=j+1; i <m;i++){
				if(Math.abs(LUcolj[i])>Math.abs(LUcolj[p])){
					p =i;
				}
			}
			if(p!=j){
				for(int k=0; k<n;k++){
					double temp = LU[p][k];
					LU[p][k]=LU[j][k];
					LU[j][k]=temp;
					
				}
				int k=pivitrow[p]; 
				pivitrow[p]=pivitrow[j];
				pivitrow[j]=k;
				pivit = -pivit;
				
				
			}
			
			if(j < m &LU[j][j] !=0.0){
				for(int i =j+1; i<m;i++){
					LU[i][j]= LU[i][j]/LU[j][j];
				}
			}
			
			
		}
		Matrix A = new Matrix(LU);
		return A;
		
		
	}
	 
	
	  public Matrix getL () {
	      Matrix Lower = new Matrix(m,n);
	      double[][] L = Lower.getArray();
	      for (int i = 0; i < m; i++) {
	         for (int j = 0; j < n; j++) {
	            if (i > j) {
	               L[i][j] = LU[i][j];
	            } else if (i == j) {
	               L[i][j] = 1.0;
	            } else {
	               L[i][j] = 0.0;
	            }
	         }
	      }
	      return Lower;
	   }
	
	  /**
	   * 
	   * @return Upper
	   */
	  public Matrix getU () {
	      Matrix Upper = new Matrix(n,n);
	      double[][] U = Upper.getArray();
	      for (int i = 0; i < n; i++) {
	         for (int j = 0; j < n; j++) {
	            if (i <= j) {
	               U[i][j] = LU[i][j];
	            } else {
	               U[i][j] = 0.0;
	            }
	         }
	      }
	      return Upper;
	   }
	  
	  public Matrix solve (Matrix B) {
	      if (B.getRowDimension() != m) {
	         throw new IllegalArgumentException("Matrix row dimensions must agree.");
	      }
	      for(int j = 0;j < n; j++){
	    	  if( LU[j][j] == 0 ){
	    		System.out.println("No solution");
	    		return null;
	    	  }
	    	  
	      }
     
	      // Copy right hand side with pivoting
	      int nx = B.getColumnDimension();
	      Matrix Result = B.getMatrix(pivitrow,0,nx-1);
	      double[][] X = Result.getArray();

	      // Solve L*Y = B
	      for (int k = 0; k < n; k++) {
	         for (int i = k+1; i < n; i++) {
	            for (int j = 0; j < nx; j++) {
	               X[i][j] = X[i][j]- X[k][j]*LU[i][k];
	            }
	         }
	      }
	      // Solve U*X = Y;
	      for (int k = n-1; k >= 0; k--) {
	         for (int j = 0; j < nx; j++) {
	            X[k][j] /= LU[k][k];
	         }
	         for (int i = 0; i < k; i++) {
	            for (int j = 0; j < nx; j++) {
	               X[i][j] =X[i][j]- X[k][j]*LU[i][k];
	            }
	         }
	      }
	      return Result;
	   }
	  
	
	/**
	 * @param args
	 */
	public static void main(String[] args) {
		// TODO Auto-generated method stub
//		double[][] array ={{1,2,3},{14,5,6},{7,8,9}};
//		double [][] result ={{12},{16},{28}};
//		System.out.println(Arrays.deepToString(array));
//		System.out.println(Arrays.deepToString(result));
		
		Readarray read = new Readarray();
		List<Matrix> list = read.readFile("linearsolver_hyin_input_002.txt");
		
		
//		Matrix rs = new Matrix(result);
//		Matrix matrix1 = new Matrix(array);
//		Matrix prove = matrix1.solve(rs);
//		System.out.println("prove"+ Arrays.deepToString(prove.getArray()));
		Matrix A=list.get(0);
		Matrix result=list.get(1);
		
		LuDecompostion test = new LuDecompostion();
		test.guassianElilmination(A);
		
		Matrix L =test.getL();
		
		double[][] Lowerarray = L.getArray();
		DecimalFormat f = new DecimalFormat("#######################################.00");
		for(int i=0; i < Lowerarray.length;i++){
			for(int j =0; j < Lowerarray.length; j++){
				
				Lowerarray[i][j] = Double.parseDouble(f.format(Lowerarray[i][j]));
			}
		}
		
		System.out.println("Lower Matrix(L) = ");
		for(int i =0; i< Lowerarray.length;i++){
			System.out.println(Arrays.toString(Lowerarray[i]));
		}
		
		
		Matrix U =test.getU();
		double[][] Upperrarray = U.getArray();
	
		for(int i=0; i < Upperrarray.length;i++){
			for(int j =0; j < Upperrarray.length; j++){
				
				Upperrarray[i][j] = Double.parseDouble(f.format(Upperrarray[i][j]));
			}
		}
		System.out.println();
		System.out.println("Upper Matrix(U) = ");
		for(int i =0; i< Upperrarray.length;i++){
			System.out.println(Arrays.toString(Upperrarray[i]));
		}
		
		System.out.println();
//		System.out.println("L"+ Arrays.deepToString(L.getArray()));
//		System.out.println("U"+Arrays.deepToString(U.getArray()));
		Matrix C = test.solve(result);
		DecimalFormat f2 = new DecimalFormat("#########################################.000");
		
		if(C!=null){
			double[][] rs = C.getArray();
			for(int i=0; i < rs.length;i++){
				for(int j =0; j < 1; j++){
					
					rs[i][j] = Double.parseDouble(f2.format(rs[i][j]));
				}
			}
			System.out.println("Solution: "+Arrays.deepToString(rs));
			
			
		}
		


		


		
	}

}
